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A description of photoassociation by CW laser is formulated in the framework of grid meth- 
od ■ ods. The Hamiltonian describing one or several bound states coupled to a multiple of continuum 
^ , manifolds via a radiative field is written in the energy representation and diagonalized. The gen- 
erality of the treatment allows to compute accurately and efficiently physical properties such as 
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non-linear high-intensity energy shifts, line shapes, and photoassociation rates both for isolated 
and non-isolated resonances. Application is given to sodium photoassociation in the experimental 
conditions of Mc Kenzie et al [Phys. Rev. Lett. 88, 090403 (2002)]. Inverted region for the 

> 

C*~> ' dependency of the rate vs. the intensity and non-symmetric lineshapes were predicted to occur 

above the saturation limit. Comparison with the model of Bohn and Julienne [Phys. Rev. A 60, 
O ; 414 (1999)] i s discussed. 
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I. INTRODUCTION 



Many experimental groups are presently working on the formation of ultracold molecules 
and molecular quantum gases. A very efficient scheme is the photoassociation reaction, 
where two colliding atoms absorb a photon to make a molecule moreover ultracold 
photoassociation spectroscopy [2j is capable of providing deep insight into the properties 
of long range molecules, cold atomic collisions, or pair correlations in a condensate. The- 
oretical two-body scattering models have been developed to compute the photoassociation 
rates , in the framework of a perturbative treatment, adapted to situations 

where low-intensity CW lasers are employed and many-body effects are negligible. A linear 
variation of the photoassociation rate as a function of intensity is then predicted. In the 
work of Bohn and Julienne 6|, which develops a semi-analytic theory yielding laser-induced 
energy shifts and 1 ine broadening, the departure from this linear behavior is described, the 
saturation limit at larger laser intensities being evaluated. In a condensate, the two-body 
model should fail and Javanainen and Mackie [§] have developed a many-body theory 
which predicts a saturation limit occurring at lower intensities and attributed to rogue 
photodissociation. The search for saturation effects has stimulated several experiments. 
Measurement, by Mc Kenzie et al joj], of the photoassociation rate in a sodium condensate 
as a function of laser intensity has shown a linear behavior in the intensity range considered, 
in agreement with the two-body Bohn and Julienne estimation, and in contrast with the 
prediction of Javanainen and Mackie. For a quantum degenerate lithium atomic gas, 
measurement [kJ of the intensity dependence of the photoassociation rate has verified both 
the large value of this rate (at a temperature T=600 nK) and the two-body saturation 
effect at intensities / > 30 W/cm 2 . The agreement with theory is considered as reasonable 
within the present experimental accuracy: the higher value of the measured saturation 
limit, and the oscillations in the experimental rate in the 30-80 W/cm 2 intensity range 
should be confirmed by measurements with a smaller error bar. For a non-degenerate gas of 
cesium atoms, at a temperature of 40 p,K, within an optical dipole trap, saturation limit is 
reached for intensities above 100 W/cm 2 |ll|], and the experimental values then lay around 
the theoretical curve : it is not clear whether such oscillations are due to uncertainties or 
should be attributed to a physical effect. 
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Recently Naidon and Masnou-Seeuws [12j have developed a time-dependent many-body 
theory of photoassociation in an atomic Bose-Einstein condensate, showing that several 
physical effects depend on pair correlations. They show that the two existing models for the 



description of such correlations 



141 ]. differ in their predictions of the photoassociation 



line shapes at large intensities. Precise measurements of the variation of the line shapes as 
a function of cw laser intensity should therefore give insight into the optimal description of 
many-body effects. 

The aim of the present paper is to anticipate such measurements by revisiting the 
theoretical treatment of photoassociation line shapes, taking advantage of the efficiency 

nnnn 

of grid methods [15J, Ha, H7|, [18( for the numerical description of cold collisions and long 
range molecules, particularly in cases where several channels are involved. Actually, 
;he later methods have been widely used for the computation of photoassociation rates 
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y, 21, £ 



22], since they can determine very precise energies and wave functions, 
both for bound and for continuum levels. It has been demonstrated that due to their 
global character, mapped grid methods are ideally suited to precise computation of the 
perturbations occurring in molecular spectra, such as Rb2 23, 24, 25] or K 2 j26] singlet- 
triplet mixing. In such cases, the perturbation of one level of a given series should be 
attributed to a large number of bound or continuum levels of the other series, and local 
perturbation treatments are poorly adapted. However, up to know, attempts to compute 



the photoassociation line shapes via grid methods have not been very efficient |27l . 128 ] 
and the interpretation of precise measurements is routinely performed with Numerov 
type approaches 29j. The present paper is proposing an alternative way of implement- 



ing grid methods for the calculation of photoassociation line shapes, energy shifts, and rates. 



The paper is organized as follows: The description of the present model is reported 
in Section II-IV. The link with previous models, such as the Bohn and Julienne model is 
discussed in Section V, where the various assumptions and approximations are underlined. 
Application to the cases of Sodium photoassociation in the experimental conditions of Ref. 
|9j is presented in Section VI. Section VII will summaries and discuss the outlook of the 
method. 
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II. THE PRESENT MODEL - PRELIMINARY REMARKS 



The physical process of photoassociation (PA) can be divided into three steps. Initially 
(step I), a thermal ensemble of cold atoms is in the ground electronic state; the probability 
to find a colliding pair with relative motion at energy E and angular momentum I is given 
by: 

M E) = '^"^W, (1) 

where T is the temperature, and is the Boltzmann constant. The partition func- 
tion Z eq was shown in [30j to be well approximated by the classical expression 
Z eq = QtV = V(2nfxkBT) 3 / 2 /h? for a gas of non-interacting pairs in a volume V. 
Further discussion regarding Z eq could be found in [scj. 

The radial wave function which corresponds to the pair is represented with a continuum 
eigenfunction \i^%i{R)} of the single-channel Hamiltonian H sc 

h h2 d2 i una i vim (2) 

Here, /i is the reduced mass of the pair, and V (R) is the Born-Oppenheimer potential curve 
for the relevant, here the ground, electronic state. 

At the second step (step II), the intensity of a CW laser with frequency u is ramped 
within a time scale that can be varied from hundreds of ns to several /is. Along this 
preparation period, the continuum states are coupled to one or more bound vibrational 
levels v, J that belong to one or more excited electronic states e. Each of the levels in the 
excited state, denoted by |0^j(i?)), may decay spontaneously with a natural free width 
7° j, which is typically at the order of ns. The field that is applied on this stage varies with 
time, and so are the states of the total system. 

Finally, during the stable period of the experiment (step III), the laser reaches its peak 
constant intensity /. For the trap-loss or multiphoton detection scheme which will be de- 
scribed here, the duration of this period is on the order of ms. 
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Each of the energy levels of the system at this step can now be written as: 

e,v,J I v',l 

where we include the possibility of populating bound levels \i>^ t ) of the ground state. In 
the present work, we shall consider only this third step, which of course is influenced by the 
preparation step. Since we shall not treat explicitly the preparation step, and possible tran- 
sitory effects, a rapid switching on of the photoassociation laser will be assumed, justifying 
sudden approximation. The density of probability for a transition from initial state \iI>ei) 
to a final state \^fj(E)) is given by: 



At section IV it will be shown that a knowledge of P(E,i)^j is sufficient to extract the 
experimental parameters of the PA process. An efficient method to compute the initial and 
the final states, \^ g E ,i) an d \^j(E)) : is the subject of the next section. 

III. NUMERICAL DESCRIPTION OF THE FIELD-MOLECULE SYSTEM 

The initial state is a thermal ensemble composed of atom pairs, with several angular 
momenta /. Each pair is coupled to several J partial waves components on the excited 
state. For simplicity, the equations will be written hereby for one value I andJ of the 
angular momentum. To obtain P(E)^j one has to compute the initial and final eigenstates 
of the Hamiltonian without and with the light interaction. In the global grid method used 
here, this is achieved by representation of the Hamiltonian in a certain basis, and then a 
diagonalization. The diagonalization step is the most expensive computationally, and it 
scales as N 3 where N is the number of basis function. A selection of an appropriate basis 
prior to the diagonalization is hence crucial for the efficiency of the method. 
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A. The initial states: ^-representation of the unperturbed molecule in a box 



Using numerical grid methods, a numerically complete set of of eigenvalues and eigen- 
vectors of the time-independent Schrodinger Hamiltonian H^, is obtained for the isolated 



molecule, for the ground and excited potentials. 



Fourier grid method (see for example the references 15|, lla, Ll7|, l31|] ) . Within the procedure, 



his is achieved by using the mapped 



a unitary transformation XJ mapp is performed on the Hamiltonian eq. (jSJ), which is written 



\Jl napp H.fJj map p, is represented on 



in equally space grid R. The resulted Hamiltonian Hf c 
a grid with non-uniform distribution of points in the position space x, to take advantage of 
the phase-space occupancy of the problem. The potential energy operator is written simply 
as V(x), and the kinetic energy operator reads: 

-2 



7T 



J V2f)J lDJ-1/2 



(5) 



L x is the length of the grid, D is the first derivative operator matrix, given by (assuming 
even number of grid points): 



D 



'•j 



d 

dx 



0,i = j 

(-ir j 

sin[(i— j)ir/N] 



(6) 



and J (x) = is the (diagonal) Jacobian operator for the transformation. 

A sufficient representation for H^. by an equally spaced grid demands typically N ex 10 4 
basis functions for the example considered below in section IVll As was shown in 17|], the 
use of a mapped grid leads to a reduction of about an order of magnitude in the basis size, 
to several thousands. 



A diagonalization of Hf c gives the initial, interaction- free, wave functions and their 
energies for both electronic states \ip 9 E {R)) and \(j>l j(R))- The eigenvalues with E < 0, 
correspond to bound, discrete states and the ones with E > refer to box-normalized 
quasi-continuum states. The continuum under this picture is therefore tied to the box size 
from which is was deduced. 



In order to take account of the spontaneous emission in the excited channel e, an imagi- 
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nary potential is introduced: [32] 

V d (R) 



2 

3c 



\n eg {R)\ 2 u{Rf 

h 



(?) 



where fi eg (R) is the transition dipole moment, and u (R)2tt (V e (R) — V g (R)) /h. 

A diagonalization of H^ c + results a complex eigen-energies. The imaginary part of 
the eigenvalues, 7° j = Im (-Eyj) is interpreted as the lifetime of the level (R))- 



B. The final states 

1. The field-dressed molecule 

A CW laser with constant intensity /, detuned by A a relative to the atomic resonance 
line, introduces a radiative coupling between the ground state and the excited state. After 
making the rotating wave approximation, the full dressed Hamiltonian under the interaction 
can be written as: 

The coupling between the ground and the excited state is introduced by the Rabi coupling 
operator through the transition dipole moment: £l(R) = —AjX eg {R)^J where fi eg is 
the spatially dependent transition dipole moment between the two electronic states, and 
A is an angular factor depending upon the polarization of the laser and the orientation 
of the molecular axis. In the R representation Q(R) is a local, i.e. diagonal, operator. A 
diagonalization of gives the final states \^j(E)). 



Hf c Cl(R) 
Cl(R) H^ c + A Q 



(8) 



2. Energy eigenstates representation 

For the final states, one has to diagonalize the Hamiltonian describing the interaction of 
the molecule with the radiation field routinely for many experimental parameters. A direct 
diagonalization of Hr even for the mapped grid representation, is inefficient, mainly due to 
mainly due to the redundancy in the representation of the basis for both potentials. An 
improved representation for the coupled system is hence crucial. This can be obtained by 
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writing in the energy eigenstates of the single channels Hamiltonians H sc . The unitary 
transformation U^_ rep is performed on to give 



H 



E-rep — U^_ rep H#U£_ rep , 



(9) 



the Hamiltonian in the energy representation. The eigenvectors that were obtained as the 
solutions for the single channel problem serve as the transformation matrices between the 
two representations, so that: 



U 



E—rep 



fl 9 

^ E-rep 




TP 

E-rep 



(10) 



where U|(f are the matrices which contain the complete set of eigenvectors for the 
ground and excited states, \i/j 9 e (R)) and \4>% j{R)), respectively. 



The coupled channel Hamiltonian in the energy states picture becomes: 



• for the diagonal part in the excited channel subspace 

(R) \ H \4>l (R)) = E° v e + A a - * 7 ° e , E° v e < 0, (11) 

where A a is the detuning from the atomic resonance, and 7° e the natural width for 
the vibrational level e, v. 

• for the diagonal matrix elements in the ground channel subspace: 

(r n (R)\m 9 n(R)) = E° n 9 - (12) 

For clarity, explicit / dependency is ommited. We remark here again that E° n 9 include 
both the bound and the quasi-continuum parts of the spectrum. E^ 9 > is therefore 
discrete, for a given box normalization. 

• for the non-diagonal part of H£_ rep : 

hn e g % = (<pUR)\m 9 n (R)) (is) 
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We can employ the R-centroid approximation for fi eg (R) and write: 



(14) 



where now Tg )V n denotes the Frank-Condon factor between vibrational levels from the 



two different electronic states. 

The advantage of the energy representation is the possibility to pick the energy levels 
that participate in the process and thus reduce the required basis set. We shall consider 
formally the case of an isolated resonance where the ground state continuum \tp^ (R)) 
is coupled only to single bound level in the excited state \(j)%{R)). The detuning from 
the atomic resonance A a could be replaced by A, the detuning from the excited bound 
level. A generalization to several bound excited states, or even a whole continuum, is 
straightforward. From this point on, we reduce all the indices for the excited state. 

By moving to the energy representation we reduce the representation for the coupled 
channels problem by about one half, to the order of N « 1000 in the example discussed 
below in section [VTl 

3. Getting out of the box - from E to E' representation. 

Note that so far, the representation of the continuum is done using a box unity- 
normalized wave functions. This box-dependent basis for the continuum is arbitrary, and 
therefore it could be reduced. 

In order to do so, one had to defined a box-independent Frank-Condon factor. Box- 
independent quantities for the continuum are obtained by moving from the box-normalized 
to the energy-normalized wave functions. The relation between the two is given by: 



where the inverse of the level density, ^P-, is given simply, for a large enough box, by 
the energy difference between two adjacent quasi-continuum levels, given under the box- 
normalization constraint. An energy normalized Frank-Condon profile: 




(15) 




(16) 



9 



could now be obtained. T e {E) is a box-independent, numerically continuous function of 
the scattering energy. In fact, T e {E) is a physical measure for the energy dependence of 
the interaction and it thus the function that one would like to sample for obtaining the 
minimal energy representation. The alternative basis has to capture the functional energy 
dependence of the Frank-Condon profile, using the smallest basis. This could be done in 



principle with any sampling procedure 33[. The technical details regarding the numerical 



procedure we employed in this work is described at the appendix. 

Assuming such a minimal representation basis was found, we will denote it by \g,n') 
for the ground state continuum, for n' > % being the number of bound states. The 
expression for the quasi- continuum diagonal part of the full Hamiltonian (eq JT2j) reads now: 

(9, n'\ H \g, n') = E%j, n' > n B (17) 

and all the other diagonal parts remained unchanged. The nondiagonal part in the alterna- 
tive basis is defined through the energy normalized Frank-Condon profile as: 

/ arp 1 

^ = = T\E). (18) 



A diagonalization of the Hamiltonian (eqs. [T71fT8] ) gives the set of eigenstates for the 
coupled Hamiltonian, \^/j(E)), which include components on both channels, with possible 
contributions from all the bound or free levels. The energy eigenvalues which are obtained 
are again complex and 7 J = lm(Ej) is the width of the level \^j(E)). The transition 
probability from any initially populated state \g,n') to any final state \^j(E)), is simply: 

*V)^- = K*i(£) \9,n')\ 2 (19) 

It is important to notice here that because \^fj) is an eigenstate of the coupled Hamiltonian, 
a pair that reaches this state at the beginning of the final step of the experiment remains 
in the same state without any significant dynamics beside an exponential decay with the 
lifetime 7 J . 



10 



Representation 


N 9 




R - Equally spaced position grid 


10 4 


10 4 


x - Mapped position grid 


10 3 


10 3 


E - Energy eigenstates - boxed 


10 3 


< 10 


E' - Energy eigenstates - non boxed 


~ 10 2 


< 10 



TABLE I: Summary of the different representations used in this paper. N g and N e are the typical 
basis set sizes for the ground and the excited electronic state, respectively, adapted to the example 
of section IVII 

Table 1 summarizes our procedure so far with the different representations and their 
typical necessary basis sizes. We started from the Hamiltonian on the equally spaced 
grid representation. A minimal representation of the problem with such a grid demands 
a basis of about N R = 10 4 function. We then used the phase space occupation of the 
problem to move to a mapped grid, and reduce the basis set in about an order on 
magnitude to N x = 10 3 , for each of the electronic states. The transformation to the 
energy representation reduces the excited state manifold. The reduced E' representation 
minimizes the basis for the ground state, so that finally we could solve the whole coupled 
channels problem with a basis at the order of Ne> ~ 200, or even smaller. At this sizes of 
basis the method become efficient, even compared to the semianalytic perturbative methods. 

In the next section the deduction of the various PA experimental parameters from P( Si7l q_ + j 
(and 7 J ) is presented. 

IV. EXTRACTION OF THE ENERGY SHIFTS, TRAP-LOSS RATES, AND 
LINESHAPES. 

With the formal definitions of the previous subsection we can now move to define the 
various physical quantities that are being measured on the system. The probability to find 
an atom pair at the beginning of the stable period with energy Ej (which we will denote as 
zero time) is given by: 

N 2 (0) 

p(Ej, 0) = N pmr (0) P( 9 ,n>)->jPB(E n ,) = at °™ K ' P( 9 ,n')^PB(E n >) (20) 

n' n' 
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where N atom (0) and N pair (0) are the number of atoms and atom pairs at time 0, and pb is 
defined above (see eq. ([I])). As was explained at the last paragraph the only dynamics at 
the stable period are the exponential decay of the j-th level. At later times thus: 

p(Ej,t) = exp(- 7j t/%(^,0) (21) 

so that the total probability density for the whole ensemble is: 

N 2 (0) 

Ptotit) = $>(£„ t) = -^fU. eM-imPM^PBiEn') (22) 

3 n',j 

The observed reaction rate v(t) at any given time is the change of Ptot{t) versus time: 
v(t) = -^f^ = ^^^IjeM-ljt/m^PBiE^) = ( 7 (t)> (23) 

n',j 

The rate could be thus interpreted as the thermally averaged lifetime. On the other hand, 
the experimental definition of the diatomic reaction rate for a sample containing N atoms is: 

N 2 (t) 

v{t) = -K at °™ yi) . (24) 

Now, a comparison between the two expressions yields the relation between the calculated 
and experimental rate constants, which holds for any given time. Specifically for t=0 one 
gets: 

K{I ,Q = vmi = m ,25) 

In a fashion similar to the experimental, a calculation of if as a function of the detuning 
yields the lineshape. The peak value of K as a function of A, and the value of the detuning 
at the peak are termed as the maximum reaction rate K max and the energy shift E\ for a 
given intensity and temperature, respectively. 



Before moving forward, we remark that the model is non perturbative with respect to the 
laser intensity. It is also well oriented to dynamical, time dependent-type calculations. Fur- 
thermore, a generalization of the method to multiple coupled continuum manifolds, several 
bound states, or to include hyperfine structure, is strightforward. Its drawback is mainly 
being a little less efficient computationally, mainly due to the expansive diagonalization 
step. Nevertheless, this could be overcomed, especially for calculating lineshapes and satu- 
ration effects, as we will show at the application of section VI. But, before doing so, we will 
present briefly in the next subsection the necessary expressions of the method of Bohn and 
Julienne^ (B & J) to which our expression (125]) will be compared to. 
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V. ADJUSTMENT OF BOHN AND JULIENNE MODEL TO THE GRID METH- 
ODS FRAMEWORK 



In this subsection we present the method of B & J [6J, and adjust it to take advantage 
of the available grid methods. We will limit the treatment to the simple case of a single 
isolated resonance embedded in a continuum. The scattering-theory-based method defines 
a scattering amplitude Se for a pair of atoms in the center of mass framework. The square 
of Se is referred to as the probability to initiate in a unbound state with kinetic energy E 
and be lost at infinite time, due to a light-induced coupling to a bound molecular state. 

According to Fanof" 



34j . \Se\ 2 is given by: 

\q 12 _ ,.)<;, 



here A is the detuning of the light frequency of the resonance bound state's energy , and 7 is 
the natural linewidth of the bound level. E\ and Qe are the energy shift and the stimulated 
absorption/emission rate, respectively. With the help of second order perturbation theory 
it was shown in (3^] that: 

Ei = V.J dE'^ (27) 

where V. stands for 'the principle part of, and Qe is the dipole moment coupling that was 
defined above (see eq. ffTB"]) ). 

The rate of loss from the system is given by performing a velocity thermal averaging of 
\Se\ 2 , for all the initial scattering energies, and with respect to all the relevant partial waves, 
to obtain the observed rate loss: 

k = "£t f PB ^)\ s ^ dE (28) 

1 

where the integral over the energy indicates for a thermal averaging over the Boltzmannian 
distribution p B , defined above (see eqJTJ. The superscript for S denotes the rotational 
quantum number for the partial wave I, and v and k = ^/{2fiE)/h are the relative velocity 
and the wavenumber of the atom-pair. 

It is important to comment here, that in order to calculate the energy dependent Qe for 
eqs. (26-28) the energy dependency of the continuum has to be calculated [see eq. (0Q)]. 



The common use of local solvers like the Numerov propagator 33j for that purpose made 



Fano's formalism impractical computationally. An energy independent Qe was assumed for 
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demonstrative purposes in Ref . 34j . Obviously, as one can see from figure 3 this is far from 
being adequate to describe molecular systems in any other circumstances. 

Using the ideas of Du et al. 35| Bohn and Juliene combined the couple of regular and 
irregular solutions of eq.(T2]) to write the energy shift E x as: 



oo R 

E x = -1tx J r Q E (R)dR J^e (R')dR' 





(29) 



where the superscripts denote the regular and irregular part of the continuum wavefunction 
for E — > 0, and r l % VLE(R) is understood to be the integrand of the integral of eq. (1131) . As 



29 



noted in 
kernel G = (^H — E 
Fano's formalism. 



351 ]. this alternative expression could be viewed as a way to write the Green 



in the position representation instead of the energy space used in 



To reduce even further the effort in the computation of Qe, the assumption of a slowly 
varying continuum with respect to E is employed, i.e., lim^o «C Pb(E — > 0). Under 
these conditions, the rate is assumed to be influenced by the on-resonance light interaction. 

The continuum is represented now by only a single continuum state and the thermal 
averaging in eq.(l28l) can be avoided. As shown in[36|], this leads to: 

2/ifc( 7 + ft(£ R ))2 1 ; 

where now the coupling fl E (defined above in eq. ffTBl ) is taken only on the resonance energy 

Er. 

An analytic alternative way for calculating the energy dependence of Q was suggested 
in Ref. |37|. In the latter work, strong energy dependence, was found in some cases, for 
instance 85 Rb. 

In this paper, we adopt the advantages of the grid methods in the E-representation, and 
calculate the explicit dependency of Q(E) directly from the continuous energy normalized 
Frank-Condon profile with any desired accuracy for any energy range. This allows us to 
calculate Ei directly using Fano's formalism, without a need to have the irregular solution 
of the single channel Schrodinger equation. As noted in 29( the two expressions give the 
same results. 
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At the next section a comparison between the various expressions will be demonstrated 
by a concrete application. 



VI. APPLICATION 



The example chosen for numerical application will be the on-resonance photoassociation 
of the v=135 , J—l level of Na2 A 1 S+, in the experimental conditions of the NIST experi- 
ment The initial state is a pair of colliding atoms with I = or 2 (see figured]). The 



calculations have been performed using the same potentials as in Ref. [121] . The position of 
the repulsive wall in the ground state X X S+ has been slightly modified, so that the binding 
energy of the last levels is E(v = 65) = 0.013cm _1 for Z=0 and E(v = 64)=0.25cm _1 
for 1—2. The experimental values for these levels are E(v = 65) = 0.013cm -1 and 
E(v = 64) = 0.25cm _1 , respectively. In the excited curve, the level t>=135 is bound by 
49.23 cm -1 , the two neighboring levels being v—134 at -52.95 cm -1 and v—136 at -45.76 
cm -1 . Note that in experiment the binding energy of the resonant level is 43 cm -1 . The 
natural width is set to 7 = 2 ttx 18.36 MHz to match the experiment. 

Figure [2] presents the coupling profile, i.e., the Frank Condon factors multiplied by the 
Rabi frequency (see eq. II4"|) . for various bound levels n of the ground state coupled to one 
level (v=135) of the excited state at intensity of 1 W/cm 2 . In the main frame n varies from 
to 64 and oscillations in the coupling term are visible as the Condon point is moving to 
large internuclear distances . The FC factors are identical for both partial waves when 
n < 64 . The n = 65 level is bound only for the S wave. 

At the upper and the lower panels of figure [3] the coupling profile for the continuum 
part of the spectrum is shown, for both S and D partial wave, respectively. The left and 
the right panels show the same profile in logarithmic and linear scale. The last bound level 
n = 65 for the S wave, is shifted for the D wave and appears as a shape resonance with 
a well defined width in the D continuum. For the continuum part of the spectrum, the 
Condon point is energy independent. The structure of the coupling profile of the continuum 
is a signature of the energy dependence of the continuum wave functions [30]. The rotational 
barrier on the D potential surface makes the D profile to vanish more strongly near the 
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Ground Electronic State 



FIG. 1: A schematic view of the CW laser PA experiment. Atoms pairs colliding in S (I = 0) 
and D {I = 2) partial waves on the ground electronic state, are coupled by a laser with constant 
frequency and intensity to the v = 135, P ( J = 1) bound level in the excited electronic state. The 
nearby bound levels on the excited state, as well as few bound levels within the ground electronic 
state could be coupled as well. The natural lifeltfime of the bound state 7 is 2ttx 18.36 MHz. 
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FIG. 2: Variation of the coupling term between the level v= 135 of Na2 A E| and the vibrational 
levels of the Na2 ground state vs. the binding energy. The intensity of the CW laser is 1=1 
W/cm -2 . The rotational quantum number is assumed to be either or 2. The first 64 levels are 
shown. Inset: the same for the upper levels for a rotational quantum number 1= (full rectangles, 
n=64 and n^=65) and 1=2 (striped rectangle, n = ra ( g=64). 

threshold than does the S profile. 

In the next subsections the results of grid method will be presented in three steps: First, 
the light induced energy shift will be calculated and discussed. Then, the results for the rest 
of the PA observables will be presented in the low intensity regime, where the perturbative 
method should be adequate. Finally, the results for in the high intensity regime, where the 
appearance of saturation effect are explored. 
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A. Light Induced Energy Shifts 

Using Fano's perturbative expression for the energy shift (eq f27|) and the right panels 
of figure [3] it is understood that the energy scale which is needed for the calculation of the 
energy shift is larger in many order of magnitude then the thermal energy, (see the appendix 
for a numerical example). On the other hand, as will be shown below, the lineshape is 
influenced only by the coupling profile in the vicinity of the thermally populated energy. 
This separation of scales can help us furthermore in reducing the basis set, depending upon 
which physical quantity one would like to calculate. 




FIG. 3: Left: Variation of the coupling term between the level v= 135 and the energy normalized 
S (upper) and D (lower) continuum levels of the Na2 ground state vs. the scattering energy at 
logarithmic scale. The intensity of the CW laser is 1=1 W/cm~ 2 . Right: the same as the left 
panels, now in a linear scale. 

The light induced energy shift for intensities up to 100W/cm 2 , are presented in figure HI 
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The energy shift was calculated using the S continuum part of the ground state manifold. 
The green line is the energy shift that was calculated using the modified Bohn and Julienne 
perturbative method and is temperature independent by definition. In red, blue and black, 
we present the energy shift using our method for T = 20fiK, 2fiK and 200nK, respectively. 
All the curves are linear as expected. A difference exists between the energy shifts, but it is 
too small to be measurable. The difference in the slopes between the present calculations 
and the perturbative method is also within the experimental accuracy. An important result 
we observed in our calculation is that the perturbative linearity of the energy shift with 
respect to the intensity was conserved for any feasible intensity that was checked (up to 
100kW/cm 2 and higher). Another feature of the perturbative limit is the additivity of the 
contributions from different partial waves with respect to the bound and the continuum. 
Figure [5] presents the energy shift for the full interaction, i.e., bound and continuum part 
of the spectrum for both S and D partial waves. No deviations from the perturbative 
treatment were found with respect to the additivity. 

B. Low Intensity regime: the perturbative limit 

Saturation effects for the case of photoassociation in a sodium condensate were not 
observed up to 1.2kW/cm 2 . A typical lineshape for an intensity of 100W/cm 2 , well below 
the saturation limit, is presented in Figure [6j The lines were all centered to zero detuning. 
The red line presents the rate versus the detuning for our method for T=500 nK. In black 
line the correspondent result for the Bohn and Julienne method with the same temperature 
is shown. The two lines show the same width but differ at the peak height by less then 
10%, small comparable to the experimental accuacy. Moreover, as can be seen from the 
blue line in the figure, a small change of the temperature for the perturbative method 
(to 400nX) reproduces exactly the same lineshape as in our method. Figure [7] presents 
with the same colors of figure [6] the results for the maximum rate lineshape K max as 
a function of the intensity. All the curves are linear, as expected for the perturbative 
regime. Again, the agreement between the two methods is within the experimental accuracy. 

We have to remark here that, to our best knowledge, an extensive study of the 
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FIG. 4: Energy shift as a function of the intensity. The values for the perturbative method are in 
green, irrespective on the temperature. In red, blue (barely visible) and black: non perturbative 
method, for 2fiK, 20fj,K, and 200nK, respectively. Only the interaction with the S continuum is 
considered. 

temperature dependency of photoassociation lineshapes has never been performed. Figure 
[8] presents the temperature dependency of K max for intensity of 100W/cm 2 . The red and 
blue lines present the values for the thermal averaged rates of expressions (|25f28|) . The 
difference between the methods at the most common temperature for PA experiments in 
condensates, i.e., a few hundreds of nano Kelvins to micro Kelvins, is at the order the 
experimental uncertainty. The values deviate significantly from each other at relatively 
high (T > lOfiK) and low (T < 100ni^)temperatures. The black line in figure [8] represents 
the values obtained from the perturbative method without energy averaging (see eq. [30]) . A 
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FIG. 5: The same as in figure HJ with the full interaction: S and D continuum manifolds and 
bound states. 

single state approximation for the rate seems to be inadequate, at least for the description 
of most of the temperature dependence of th e trap loss. In addition, we checked the 
results obtained applying the perturbative model (|28|) with thermal energy averaging but 
assuming a constant coupling Q(E). The cyan and green lines show the values for Q(E) 
taken at k^T and 0.5kbT, respectively. The constant coupling approximation fits perfectly 
to the fully thermally averaged model for the low temperature regime with E = 0.5kfyT. 
Two factors play a role in this situation: the coupling profile (see fig. [3]) and the Boltzmann 
thermal distribution (see eq. [1]). The first increases with E while the last decreases. From 
the results it seems that in our case, the compensation of the two factors takes place at 
E=0.5KbT. Above this temperature, the thermal width starts to play a role, and the 
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FIG. 6: A typical lineshape, i.e., rate of loss vs. the detuning, for the low intensity regime, (red) 
The non-perturbative method for 500nK. (black and blue) The perturbative method, for 500nK 
and 400nK, respectively. All the lines correspond to full energy averaged calculations, and are 
shifted to be centered at zero detuning, for comparison. 

constant coupling approximation is no longer valid. 
C. Non Perturbative regime - saturation effects 

We now move to explore the regime of high intensity. In figure [9] we present the intensity 
dependence of K max for / up to 300 kW/cm 2 . The blue, black and red lines are the 
results for our method, the thermally averaged, and the non averaged perturbative method, 
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FIG. 7: The intensity dependency of the peak loss rate values for a given lineshape (K max ) for 
the non-perturbative and perturbative methods. Colors are the same as in figure 5. All lines 
correspond to thermally averaged calculation for the assigned temperatures. 

respectively. The temperature is 200nK, typical for a condensate. All the lines show 
saturation behavior, but the saturation limit is not the same. The results for the single 
energy perturbative approximation are significantly lower then the other two methods. The 
intensities for reaching saturation are way too high to be observed experimentally, and, 
indeed, were never observed. The validity of the two-body model does not hold for such 
intensities, and therefore the significance of the results is mainly in predicting a saturation 
effect for other species where the saturation is reached at lower intensities, e.g., cesium [ll]. 
lithiumflOl] or strontium |38| . 
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FIG. 8: The peak loss rate for a given lineshape, K max , vs. the temperature for the three models. 
Red: this paper's model, Blue: the full thermally averaged perturbative, and Black: the single 
continuum energy approximation, see eq. (|30p in the text. Cyan and Green: the perturbative 
model with energy averaged but with a constant coupling Q, taken at k^T and 0.5ki,T, respectively. 



Figure [[(]] presents the same calculations as figure [9], but for T = 20/iK, a typical 
temperature for a MOT. At a MOT temperature saturation is reached at lower intensities 
for all the models. The non-perturbative method presents an inverted region in which 
the peak rate constant decrease with the light intensity. This prediction has never been 
observed, and it might be feasible experimentally for species other then Sodium. The 
saturation for the case of cesium in a MOT may give clue to this inverted behavior 11] 
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FIG. 9: Same as figure ([7]), now for high intensities and T = 200nK, a typical condensate tem- 
perature. Blue and Black and the non perturbative and the perturbative fully thermally averaged 
models. Red: a single continuum level at 200nK. 

A better understanding for the origin of the inverted behavior could be acquired by an 
inspection of the lineshapes near the threshold. Figure [11] exhibits the lineshapes computed 
with the non-perturbative method around the saturation regime. The four lines corresponds 
to the intensities of 1,2,3 and 6 kW/cm 2 in blue, red, black and green, respectively. All 
the lines are centered to zero detuning for comparison. A shoulder at the positive tail of 
the Lorentzian lineshape starts to develop as saturation is approached. Beyond saturation, 
the shoulder takes over the Lorentzian lineshape at the positive tail of the line to create an 



25 




FIG. 10: Same as figure [U for T = 20/iif, a typical temperature for a MOT. 

unsymmetrical lineshape. The power broadening at the negative detuning tail differs from 
the positive one, and the unsymmetrical form of the lineshape is more pronounced with 
the increase of the intensity. Figure [12] presents the lineshapes close to and beyond the 
saturation limit of the perturbative method, for T = 20fiK. The lines are centered to zero 
detuning. Beside the obvious power broadening, the lineshapes remain symmetrical even 
for high intensities. 

As a first order perturbative method, Bohn and Julienne model takes into account 
only single path processes. That is, the rate constant measures the population flux from 
the unbound continuum levels to the bound level in the excited state. Any flux that is 
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FIG. 11: Lineshape around the saturation limit, using the non-perturbative method. The temper- 
ature is T = 20fiK, and the blue, black, red and green lines are for intensities of 1, 2 3, and 6 
kW/cm 2 , respectively. All the lines are shifted to be centered to zero detuning. For the compar- 
ison, the line profiles are normalized to the line of lkW/cm 2 , i.e., the lines for 2,3 and 6kW/cm 2 
are divided by 2, 3 and 6. 

transfered to the bound level is lost at a rate proportional to the natural linewidth of the 
bound level. For low temperatures, the thermal energy width is narrow. The loss rate for 
the off resonance laser is determined only by the distance from the threshold. Positive 
and negative detuning are identical in that sense and the lineshape is therefore symmetric. 
For higher temperatures, the population spread is larger and the inner structure of the 
continuum start to play a role. This is the origin for the deviation from a Lorentzian 
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FIG. 12: Same as figure [TTl using the perturbative method. The temperature is T = 20[iK, and 
blue, black and red lines are for intensities of 5, 10 and 25 kW/cm 2 , respectively. All the lines are 
shifted to be centered to zero detuning. 

lineshapes which presented in figure 5 of jf| , for temperatures of lmK. This non symmetric 
form is expected to be washed out at low temperatures and high intensities, where the 
thermal wi dth is much smaller then the radiative width. 

In the non-perturbative method, all paths are being taken into account. At large intensity, 
multi-paths processes can take place. A flux could bounce from the bound level back to the 
continuum and vice verse. As can be seen in figure El the coupling profile indicates for a 
maximum of the interaction around 10~ 2 cm _1 . The shoulder of the lineshape appears at a 
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detuning equal to the maximum. Another support for this interpretation could be found in 
the energy dependence of the transition probability, P( g , n ')^j- At low intensity single path 
processes take place and the transition probability is significant only close to the threshold. 
At large intensities other parts of the profile play a role, Pf g>n n^j deviates from a symmetric 
Lorentzian and peaks also at the energy of the maximum of the interaction profile. A further 
increase of the intensity enhances even further the weight of the non-linear effects on the 
line shape, and so the lineshape at positive detuning is rolled by non linear effects. A clue 
to support this result could be seen maybe at the apparent unsymmetrical lineshapes that 
were measured for the saturation limit of cesium in a MOT (see 111]). 



VII. DISCUSSION AND SUMMARY 



In the present paper we have explored the advantage of grid methods for performing a 
non-perturbative calculation of photoassociation lineshapes, energy shifts and rates. Several 
physical properties determine the energy scales that are involved in the PA experiment and 
characterize the outcomes of the various observables. The energy shift, which is the value 
of the detuning at which a photoassociation lineshape peaks, is the parameter which is the 
less sensitive to the structure of the interaction profile and the thermal population. In order 
to calculate it accurately one has to consider an interaction energy which goes up to several 
hundreds of Kelvins and to include all the manifolds which participate in the interaction 
with the bound level. The interaction energy is independent of the thermal occupation of 
these levels. A full non-perturbative treatment for the energy shift is therefore extremely 
difficult, but fortunately, it is also unnecessary. The advantage of the grid methods in 
calculating energy shifts is minor and they could be obtained efficiently and accurately by 
perturbative methods. 



For low laser intensities, the natural linewidth of the bound level, which is usually at the 
order of hundred of MHz, or tens of milliKelvin, is the next energy scale in the hierarchy. 
The width of the lineshape and the scaling of the peak height are determined by this 
parameter at low intensities. As the intensity of the laser increases, the radiative coupling, 
or the Rabi frequency of the given scattering energy Q(E) begin to dominates, and lead to 
power broadening and saturation. 
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Notation 


Parameter 


Order of magnitude(Kelvins) 


Ei 


Energy shifts interaction energy 


10 1 - 10 3 


7 


Natural linewidth 


10 -1 - 1(T 3 


V E 


Stimulated emission rate 


Depends on A and I 


Et 


Thermal Energy 


1(T 7 - 10" 4 



TABLE II: Parameters involve in PA and their typical energy scales in Kelvins. 



The thermal energy is commonly the smallest energy scale. It ranges from about a 
hundred of nanokelvins for the coldest condensates, to a few hundreds microkelvins, for 
'hot' MOTs. For low temperatures, the narrow width of the populated energy band makes 
the coupling profile to be almost insignificant. Saturation could then be described quite well 
with the perturbative methods. As the temperature increases, the inner structure of the 
continuum is encrypted into deviations of the lineshape from a symmetric Lorentzian, and 
an inverted region of K max as the intensity increases. Some of these deviations, especially 
the ones which are visible at the low intensity limit, for relatively high temperatures, 
could be well described by the single-path perturbative models. Other features which 
result only at high intensities, could be calculated only with the help of the multi-path model. 



The non perturbative model that we presented in this paper can serve as a useful tool to 
model easily other cases where non trivial lineshapes appear. These include: 

• Photoassociation close to threshold. 

• Accidentally degenerate or almost degenerate levels that originate from different elec- 
tronic states. In both cases, interaction between several levels in the excited state need 



to be put in the 
photoassociation 



ramework of a higher order model. The case of resonant coupling 
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23l . |25| is a good example. 



Active interaction between multiple continuum manifolds, where more than a single 
continuum is populated or depopulated during the process. One possibility is when 
more then one partial wave participates in the scattering, e.g., the D resonance in 



8 5Rb (see ref. 



371]). In cases of lower barriers, higher temperatures, or in the cases 
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of processes that take place by explicitly populating the continuum, highly structured 
lineshapes might be resolved by our method. 

Finally, the dynamical approach which is the basis for the present study is an important 
building block for a future comprehensive time dependent description of the many body 
interaction in the description of PA. 
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APPENDIX 

Two purposes are the objects of this appendix: the numerical details of the sampling 
of the continuum profile T e [E) (see eq. (j!6p ). and the demonstration of the energy scale 
that is needed for the calculation of the energy shift. The sampling of the continuum profile 
T e (E) could be done in principle by a global or local numerical sampling procedure (see 
for example [331]) .The common procedures take advantage of an a-priori knowledge about 
the given function to be sampled. In this work we limited ourselves to a naive approach 
and determined the sampling of the profile a-posteriori according to the results that were 
obtained with a given basis. The criterion we employed was the numerical convergence of 
the integral over a given moment of the interaction profile T e [E\. 



where n is the order of the moment, Ej is the discrete sampling energy points, and A(Ej) 
is the energy difference between two adjacent points. The upper bound for the profile E cut 
is discussed below. The procedure initiates with a small set of points, and a new point is 
added or rejected between any two old points according to its influence on the value of X. 
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The quality of the procedure is determined by the tolerance parameter r = \z new /I old — 1 . 
The number of basis functions for the sampling the S continuum, as a function of the 
tolerance parameter is presented in figure Al. The various lines present different moments. 




Tolerence 



FIG. 13: Number of basis function vs. the tolerance parameter r. The numbers near the lines 
stand for different moments of the interaction profile J- e {E). The box normalized from which the 
original profile is deduced contains 920 basis functions. 

As can be understood, the convergence of the sampling becomes more efficient for higher 
moments, which emphasize more the features in the profile. Too high moments, however, 
smear completely the functional structure and so an optimum moment have to be found. 
We have to comment here that in order to increase the resolution of the basis set with 
respect to the thermally populated energy near the threshold (for E < NksT, N is usually 
10), we inserted a small energy grid (typically 50-100 basis functions) with constant energy 
step at close-to-threshold energies. 



35 



In figure A2 the resulted energy shifts are plotted as a function of the tolerance parameter 
for the various moment orders. The numerically exact value which was calculated using the 
perturbative method is 27.8MHz. The values for moment order above the 5 th , convergence 
to this values within less then 5% for tolerance value of 5 x 10~ 5 . The basis set that was 
taken in this paper uses r = 5 x 10~ 7 and n = 9. Larger moments than the chosen exhibit 
slower convergence and demand much smaller values of r to produce an adequate basis set. 
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Tolerence 



FIG. 14: Calculated energy shift vs. the tolerance parameter r. The intensity is WOW/ cm 2 , 
the temperature is 2[iK. Only calculation with respect to the S continuum are presented. The 
continuum is sampled up to energy of 300-fT. The numbers near the lines stand for different 
moments of the interaction profile J- e {E). 

We must remark also, that the only physical parameter which demands such a care in 
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the sampling is the energy shift. The values for the lifetime are much less influenced by 
the size of the basis, and converged even at r as small as 1 x 10~ 3 . To demonstrate the 
importance of very large energy scale that are needed to calculate the energy shift correctly 
we present in figure A3 the calculated energy shift as a function of E cut the upper bound of 
the energy grid. Convergence of the energy shift is achieved only when reaching hundreds of 
Kelvins. Note, again, that only the energy shift exhibit a sensitivity for such high energies. 
For low laser intensities the lifetime is found to converge E cut at the order of the thermal 
energy grid. Convergence near the saturation limit demands extension of the grid up to the 
first minimum of the coupling profile, i.e.,10 _2 cm _1 , or tens of milikelvins. A basis set for 
such energies contained less then 150 function and is adequate for calculating most of the 
lineshapes features, except for energy shifts. 
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Interaction Energy [K] 



FIG. 15: Calculated energy shift vs. highest continuum energy to sample. The parameters are 
identical to figure A2. The ninth moment with tolerance of 5 x 1CP 7 is taken. 
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